function res = SecondBestStep1(c1,xi,pr,pw,w,pai)
% Input  : c1  = consumption for the lowest type
%          xi  = multiplier for RC
%          pr  = parameters
%          pw  = Pareto weights
%          w   = exp(alpha)
%          pai = probablity measure
% Output : res = residual of IC for the highest type

flag = 0;
c    = zeros(pr.na,1);
y    = zeros(pr.na,1);

c(1) = c1;
phi  = pai(1)/pai(2)*(pw(1)-xi*c(1)^pr.gam);
y(1) = (xi*w(1)^(1+pr.sig)/pr.Om/(pw(1)-phi*pai(2)/pai(1)*(w(1)/w(2))^(1+pr.sig)))^(1/pr.sig);
c(2) = zero(c(1)/2,c(1)*3,eps,1e-15,@(x) SecondBestStep1IC(x,c(1),y(1),w(2),w(3),pw(2),phi,xi,pr)); 

if c(2)>c(1) % Check c2>c1
    
for i = 2:pr.na-1
    c(i)   = zero(c(i-1),c(i-1)*3,eps,1e-15,@(x) SecondBestStep1IC(x,c(i-1),y(i-1),w(i),w(i+1),pw(i),phi,xi,pr));
    y(i)   = (xi*w(i)^(1+pr.sig)/pr.Om/(pw(i)+phi-(pw(i)+phi-xi*c(i)^pr.gam)*(w(i)/w(i+1))^(1+pr.sig)))^(1/pr.sig);
    phi= pai(i)/pai(i+1)*(pw(i)+phi-xi*c(i)^pr.gam);
    if phi < 0
        res  = -c1*10;
        flag = 1;
        break
    end
end

else
   xmin = 0;
   xmax = xi*pai(1)/pai(2);
   for iter = 1:100
   flag2 = 0;
   chi   = (xmin+xmax)/2;
   phi   =  pai(1)/pai(2)*(pw(1)-xi*c(1)^pr.gam);
   y(1) = (((xi-chi*pai(2)/pai(1))*w(1)^(1+pr.sig))/pr.Om/(pw(1)-phi*pai(2)/pai(1)*(w(1)/w(2))^(1+pr.sig)))^(1/pr.sig);
   for i = 2:pr.na-1
       if abs(chi) < 1e-8
           flag2 = 1;
           ii = i;
           break
       elseif chi > 0
           c(i) = c(i-1);
           y(i) = y(i-1);
           chi   = pai(i)/pai(i+1)*(xi+chi-y(i)^pr.sig/w(i)^(1+pr.sig)*pr.Om*(pw(i)+phi-(pw(i)+phi-xi*c(i)^pr.gam)*(w(i)/w(i+1))^(1+pr.sig))); 
           phi   = pai(i)/pai(i+1)*(pw(i)+phi-xi*c(i)^pr.gam);
       else
           xmin = (xmin+xmax)/2;
           break
       end
       if  i==pr.na-1
           xmax = (xmin+xmax)/2;
       end
   end
   if flag2 == 1
       break
   end
   end
   if flag2 == 1
    for i = ii:pr.na-1
    c(i)   = zero(c(i-1),c(i-1)*3,eps,1e-15,@(x) SecondBestStep1IC(x,c(i-1),y(i-1),w(i),w(i+1),pw(i),phi,xi,pr));
    y(i)   = (xi*w(i)^(1+pr.sig)/pr.Om/(pw(i)+phi-(pw(i)+phi-xi*c(i)^pr.gam)*(w(i)/w(i+1))^(1+pr.sig)))^(1/pr.sig);
    phi     = pai(i)/pai(i+1)*(pw(i)+phi-xi*c(i)^pr.gam);
    if phi < 0
        res  = -c1*10;
        flag = 1;
        break
    end
    end
   end
end

if flag == 0
c(pr.na) =  ((pw(pr.na)+phi)/xi)^(1/pr.gam);
y(pr.na) = ((xi*w(pr.na)^(1+pr.sig))/pr.Om/(pw(pr.na)+phi))^(1/pr.sig);
if pr.gam==1
res          =   log(c(pr.na))   - (y(pr.na)  /w(pr.na))^(1+pr.sig)*pr.Om/(1+pr.sig)...
               - log(c(pr.na-1)) + (y(pr.na-1)/w(pr.na))^(1+pr.sig)*pr.Om/(1+pr.sig);
else
res          =  c(pr.na)^(1-pr.gam)/(1-pr.gam)   - (y(pr.na)  /w(pr.na))^(1+pr.sig)*pr.Om/(1+pr.sig)...
              - c(pr.na-1)^(1-pr.gam)/(1-pr.gam) + (y(pr.na-1)/w(pr.na))^(1+pr.sig)*pr.Om/(1+pr.sig);
end
res = res/10;
end

% disp([c1 res]);

save c c
save y y
